Session 5 — Statistical inference and linear regressions

Session’s plan

• Basics of Quarto documents

• From descriptive statistics to inference

• Normal distribution

• Central limit theorem

• Univariate inference

• Bivariate statistics

• Linear regressions

Learning objectives

1. Understand the basics of Quarto documents
2. Distinguish description, sampling, and inference.
3. Interpret probabilities, sampling uncertainty, confidence intervals, and p-values.
4. Choose an appropriate bivariate test.
5. Interpret linear regression coefficients and the role of control variables.
6. Understand the function of control variables
7. Know how to use R to conduct bivariate analyses and regressions

Why use Quarto?

In a normal workflow, you may have separate files:

analysis.R
figures/
results.xlsx
report.docx

This can become difficult to update.

With Quarto, the analysis and the written explanation are kept together.

What we will build today

We will start from a very small document and progressively add the pieces needed for a clean report.

analysis.qmd
     ↓ Render
analysis.html

The core idea

A Quarto document combines three ingredients.

A minimal Quarto document

---
title: "My first Quarto document"
author: "Student Name"
format: html
---

# Introduction

This is my first Quarto document.

```{r}
summary(mtcars)
```

Quarto in RStudio

In RStudio, you can create a new Quarto document with:

Then choose:

• a title,
• an author,
• an output format,
• the engine: Knitr for R.

Exercice 1

Create a Quarto document

Go in RStudio and create a qmd file with:

Give it a title, an author (you) and a .html output format

Two important buttons

Why this distinction matters

A document may work line by line but fail when rendered.

Common reasons:

• a package is loaded in the console but not in the document,
• an object was created earlier but not in the .qmd file,
• the working directory is different,
• a file path is incorrect.

The three-part structure

---
title: "My report"
format: html
---

# Text section

This is written in Markdown.

```{r}
# This is R code
mean(mtcars$mpg)
```

the YAML header

The YAML header is the part between the two --- lines.

---
title: "My report"
author: "Student Name"
date: today
format: html
---

It tells Quarto how to create the document.

YAML: the most useful fields

YAML indentation

YAML uses indentation to show that an option belongs to another option.

format:
  html:
    toc: true

Read this as:

A useful HTML header

---
title: "My report"
author: "Student Name"
date: today
format:
  html:
    toc: true
    code-fold: true
execute:
  warning: false
  message: false
---

Output formats

The same .qmd file can produce different outputs.

format: html

format: docx

format: pdf

Exercice 2

Change the YAML header of your Qurto document to include

• Today’s date
• Table of content
• Execution without warnings

Markdown: the idea

Markdown is a simple way to write formatted text.

You write plain text with a few simple symbols.

# Main title

## Section title

This is a paragraph with **bold** and *italic* text.

Headings

# Level 1 heading

## Level 2 heading

### Level 3 heading

Use headings to organise your report.

Lists

- first item
- second item
- third item

Or numbered lists:

1. first step
2. second step
3. third step

Links and images

A link:

[Quarto documentation](https://quarto.org/docs/)

An image:

![](figures/my_figure.png)

With a size:

![](figures/my_figure.png){width="60%"}

Text before code

Do not only show outputs.

Explain what the reader is about to see.

Exercise 3

In your Quarto document add:

• A first header entitled ‘Types of lists’
• A second header entitled ‘Bullet list’ with a common list as text
• A second header entitled ‘Numbered list’ with a numbered list
• A third header entitled ‘URL and images’ with one URL and one image from a URL

Adding basic R chunks

```{r}
summary(mtcars)
```

Everything inside the chunk is ordinary R code.

A setup chunk

A setup chunk is usually placed at the beginning of the document.

```{r}
#| label: setup
#| message: false
#| warning: false

library(tidyverse)
```

Chunk options

Chunk options control how the code and output appear.

Hide code, show output

```{r}
#| echo: false

summary(mtcars$mpg)
```

Useful when the reader needs the result but not the code.

Show code, do not run it

```{r}
#| eval: false

model <- lm(y ~ x, data = data)
summary(model)
```

Hide messages and warnings

```{r}
#| message: false
#| warning: false

library(tidyverse)
```

Global execution options

If you want the same option for all chunks, put it in the YAML header.

execute:
  echo: false
  warning: false
  message: false

This avoids repeating the same options in every chunk.

Figures from R code

```{r}
#| fig-cap: "Relationship between weight and fuel efficiency"

ggplot(mtcars, aes(x = wt, y = mpg)) +
  geom_point()
```

Quarto automatically includes the figure in the document.

Figure size

```{r}
#| fig-width: 7
#| fig-height: 5
#| fig-cap: "A larger figure"

ggplot(mtcars, aes(x = wt, y = mpg)) +
  geom_point()
```

Use this only when the default size does not work well.

Tables

A simple data frame will be printed as a table.

```{r}
head(mtcars)
```

For nicer tables, you can use packages such as gt, knitr, or modelsummary.

A simple summary table

```{r}
mtcars |>
  summarise(
    mean_mpg = mean(mpg),
    mean_weight = mean(wt),
    n = n()
  )
```

The important idea is simple:

Exercise 4

Create three R chunck with the code

```{r}
RSwissPos::getCAP()
```

• One not displaying the code chunck but only the output
• One displaying the code chunck but not running the code
• One displaying the code chaunk and the output

Suggested report structure

For your data driven report, a simple structure is enough.

# Introduction

# Data

# Descriptive analysis

# Model

# Interpretation

# Conclusion

Additions: Things Quarto can also do

Quarto can also handle:

• citations and bibliographies,
• equations,
• cross-references,
• websites,
• dashboards,
• revealjs slides,
• Word and PDF reports.

Equations

Inline equation:

The mean is written as $\bar{x}$.

Display equation:

$$
y_i = \beta_0 + \beta_1 x_i + \epsilon_i
$$

Citations

A citation usually requires a bibliography file.

bibliography: references.bib

Then cite in the text:

This follows previous work on reproducible research [@author2024].

Where to look things up

Useful official resources:

• Quarto guide: https://quarto.org/docs/guide/
• Markdown basics: https://quarto.org/docs/authoring/markdown-basics.html
• Quarto with R: https://quarto.org/docs/computations/r.html
• Revealjs presentations: https://quarto.org/docs/presentations/revealjs/
• Posit cheat sheets: https://posit.co/resources/cheatsheets/

What you need to remember

The minimal workflow

1. Create a .qmd file
2. Write text
3. Add R chunks
4. Render
5. Read the output
6. Fix and render again

From descriptive statistics to inference

• Descriptive statistics describe what we have.

• Inference aims to generalize from what we have to what we do not have.

Normal distribution

Normal distribution and Galton’s board

Although individual events are random, collectively we observe some patterns

If events are random, they follow a normal distribution.

What is the probability that given what we see, the events are not random?

Galton’s board

Normal distribution and certainty

• Never 100% sure of anything.

• Some threshold: 95%, 99%, 99.9%

• From the bell curve, how many SD do we need to take to be sure at this level?

Normal distribution and certainty

Central limit theorem

Even if the distribution of a variable is not normal, we reach the normal distribution if we sample the distribution and measure the distributions of means.

Central limit theorem

• Samples vary: Random samples from the same population will have slightly different means.
• Sample means follow a normal distribution: As the sample size increases, the distribution of sample means approaches a normal distribution, even if the original population is not normally distributed.
• The average remains the same: The mean of the sampling distribution equals the population mean.
• Larger samples mean greater precision: As sample size increases, sample means become more concentrated around the population mean.
• The CLT allows us to quantify sampling uncertainty, construct confidence intervals, and perform hypothesis tests.

Central limit theorem

ggplot()+geom_density(aes(qlnorm(c(1:9999)/10000, 1, .5)*10000))+
  xlab("Distribution of income")

Central limit theorem

p1 <- ggplot()+geom_density(aes(qlnorm(c(1:9999)/10000, 1, .5)*10000))+
  xlab("Distribution of income")


sample_means <- replicate(
  10000,
  mean(rlnorm(100, meanlog = 1, sdlog = 0.5) * 10000)
)

p2 <- ggplot()+
  geom_density(aes(sample_means))

ggarrange(p1, p2)

Population, sample, parameter, statistic

• Population: all units we want to learn about (e.g., eligible residents).
• Sample: observed units (e.g., 400 survey respondents).
• Parameter: unknown population quantity, often \(\mu\) or \(p\).
• Statistic: computed sample quantity, often \(\bar{x}\) or \(\hat p\).

From the CLT, we know that

\[\bar{x} = N(\mu, \frac{\sigma^2}{n})\]

Univariate inference

Descrptively, a mean of a sample is

\[\bar{x}=\frac{1}{n}\sum_{i=1}^n x_i\]

Univariate inference

With the CLT, we know that our observed mean is normally distributed on the true mean

How certain are we about our mean.

Difference between descriptive statistics and inference

• The mean of a sample is an objective value.

• We cannot really say something general about it.

• To do so, we need to quantify the uncertainty of what we observe.

• We can use the standard deviation to quantify the uncertainty and provide leels of confidence.

Bivariate analyses

• Metric X metric : Correlation

• Categorical X categorical : Chi-Square test

• Metric X Categorical : ANOVA and t-test

Bivariate analyses

• Hypothesis testing: Tests that enables us to verify/ falsify hypothesis

• Correlation –> H0: These variables are not correlated

• Chi-Squared test –> H0: The distribution of these variable are independent

• Anova/T-test –> H0: There are no difference between the groups

Correlation

Definition: Measures the strength and direction of the linear relationship between two numerical variables.

Interpretation: The correlation coefficient (r) ranges from -1 to +1, with 0 indicating no linear relationship.

Correlation in R

x <- rnorm(1000, c(1:10))
y <- rnorm(1000, c(10:1))

cor.test(x,y)

    Pearson's product-moment correlation

data:  x and y
t = -65.343, df = 998, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.9114254 -0.8878619
sample estimates:
       cor 
-0.9003012 

Chi-Square test

Definition: Tests whether two categorical variables are statistically independent.

Interpretation: A significant p-value (p < 0.05) suggests an association between the two variables.

Chi-Square test in R

pets <- sample(c("Fish", "Cat", "Dog", "Snake"), 3000, replace=T, prob=c(.08, .45, .45, .02))

gender <- sample(c("Men", "Women"), 3000, replace=T, prob=c(.4, .6))

table(pets, gender)

       gender
pets    Men Women
  Cat   534   799
  Dog   570   803
  Fish   92   148
  Snake  23    31

chisq.test(table(pets, gender))

    Pearson's Chi-squared test

data:  table(pets, gender)
X-squared = 1.2362, df = 3, p-value = 0.7443

T-test (mean difference)

Definition: Tests whether the means of two groups are statistically different.

Interpretation: A significant p-value (p < 0.05) suggests that the population means of the two groups differ.

T-test (mean difference) in R

value <- rnorm(2000, c(5,6))

group <- rep(c("A", "B"), 1000)

t.test(value ~ group)

    Welch Two Sample t-test

data:  value by group
t = -22.019, df = 1993.9, p-value < 2.2e-16
alternative hypothesis: true difference in means between group A and group B is not equal to 0
95 percent confidence interval:
 -1.0791526 -0.9026429
sample estimates:
mean in group A mean in group B 
       4.974658        5.965556 

Analysis of Variance (ANOVA)

Definition: Tests whether the means of three or more groups are equal.

Interpretation: A significant p-value (p < 0.05) suggests that at least one group mean differs from the others.

ANOVA in R

value <- rnorm(4000, c(5,6,5,4))

group <- rep(c("A", "B", "C", "D"), 1000)

model_aov <- aov(value~group)
summary(model_aov)

              Df Sum Sq Mean Sq F value Pr(>F)    
group          3   1880   626.8   647.1 <2e-16 ***
Residuals   3996   3871     1.0                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Regression models

Regressions are a statistical method used to examine the relationship between a dependent variable (Y) and one or more independent variables (X).

They estimate the direction and the strength of an association.

They enable to do statistical inference and quantify the uncertainty

They are used to estimate effect and to predict outcomes.

Regression equation

The most basic equatin for a regression is:

\[ Y_i = \beta_0 + \beta_{1}X_i + \epsilon_i\]

Visual example

Linear model

In R, we can estimate the effect of X on y with :

mod <- lm(y~x, data_test)
summary(mod)

Call:
lm(formula = y ~ x, data = data_test)

Residuals:
    Min      1Q  Median      3Q     Max 
-43.836  -9.264  -0.215   9.322  39.888 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.09084    0.84297   7.225 9.91e-13 ***
x            0.87392    0.01433  60.972  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.81 on 998 degrees of freedom
Multiple R-squared:  0.7884,    Adjusted R-squared:  0.7881 
F-statistic:  3718 on 1 and 998 DF,  p-value: < 2.2e-16

Predictions VS observations

pred <- coef(mod)[1] + coef(mod)[2]*data_test$x

ggplot()+
  geom_point(mapping = aes(x=data_test$x, y=data_test$y), size = .4)+
  geom_line(mapping = aes(x=data_test$x, y=pred), color = "red")

The error term

Each observation has an error associated

ggplot()+
  geom_point(mapping = aes(x=data_test$x, y=data_test$y), size = .4)+
  geom_line(mapping = aes(x=data_test$x, y=pred), color = "red")+
  geom_segment(aes(x=data_test$x, y=data_test$y,
                   xend = data_test$x, yend = pred), linewidth = .2, color = "green")

The general equation

\[Y_i = \beta_0 + \sum{\beta_{n}X_{n,i}}\]

Or,

\[Y_i = \beta_0 + \beta_{1}X_{1,i} + \beta_{2}X_{2,i} + ... + \beta_{n}X_{n,i}\] Regression can have as many covariate as possible.

This enables us to control for covariates.

Control variables

Why do we use control variables in regression?

• Cofounders

• Ommitted variables

• Improving precision

• Testing the robustness

A special case: Simpson paradox

We have 2 variables that are distributed as follows

Simpson paradox

We can see a positive association between the two variables.

    Pearson's product-moment correlation

data:  x and y
t = 39.226, df = 2998, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5582293 0.6055476
sample estimates:
      cor 
0.5823816 

Simpson paradox

However, a second variable x2 separates the variables into three distinct groups

Simpson paradox and control variables

• The effect observed between groups can be different than the effect observed within groups.

• To avoid this, we need to control for variables that have a potential effect on x (IV) and y (DV).

For a short summary (in French)

Elements of a regression

• Regression coefficients (positive/negative? Significant?)

• Constant (When X=0, what is the predicted value of y?)

• R Squared (Explanatory power of the model)

• Fisher test, degrees of Freedom, distribution of error

Regression coeficients

• Estimate the effect of the independent variable on the dependent variable

• When X increases, how much does Y change?

  • Direction, Magnitude and intensity

\[Y_i = \beta_0 + \beta_{1}X_{1,i} + \beta_{2}X_{2,i} + ... + \beta_{n}X_{n,i}\] In the equation, the \(\beta_n\) that multiply variables

Constant

When the matrix of covariates is equal to 0, what is the predicted value of Y?

• To have a perfect line, we need to add something to the baseline.

R Squared

Predictive quality o our model

Share of the dependent variable variance explained by the model

Between 0 and 1, where 0 mean nothing is explained and 1 mean no variance left to explain.

Distribution of errors

Indicated whether the errors are normally distributed.

Ideally, we would like our model to have as much positive and negative error (normality assumption)

We would also want the error to be distributed equally across all values of X and Y (Heterosedasticity asusmption)

Fisher test

Test whether the model is statistically sound

Whether testing this model is better than not testing it (appr. better than random guess)

Degrees of Freedom

Tells us how much degree of freedom we have in the model.

Impossible to include more variables than observations

Ideally, we would like to have models with 10 times more observation than variables.

In R

summary(mod)

Call:
lm(formula = y ~ x, data = data_test)

Residuals:
    Min      1Q  Median      3Q     Max 
-43.836  -9.264  -0.215   9.322  39.888 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.09084    0.84297   7.225 9.91e-13 ***
x            0.87392    0.01433  60.972  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.81 on 998 degrees of freedom
Multiple R-squared:  0.7884,    Adjusted R-squared:  0.7881 
F-statistic:  3718 on 1 and 998 DF,  p-value: < 2.2e-16

Presenting regressions results

• Tables

• Effect plots

• Predictions

Table with stargazer

stargazer::stargazer(mod, type = "text")

===============================================
                        Dependent variable:    
                    ---------------------------
                                 y             
-----------------------------------------------
x                            0.874***          
                              (0.014)          
                                               
Constant                     6.091***          
                              (0.843)          
                                               
-----------------------------------------------
Observations                   1,000           
R2                             0.788           
Adjusted R2                    0.788           
Residual Std. Error      13.814 (df = 998)     
F Statistic         3,717.567*** (df = 1; 998) 
===============================================
Note:               *p<0.1; **p<0.05; ***p<0.01

Effect plots

## Extract coeficients
coefs_mod <- data.frame(coef = summary(mod)[["coefficients"]][2,1],
                        sd = summary(mod)[["coefficients"]][2,2])

ggplot()+
  geom_point(aes(x=coefs_mod$coef, y="X"))+
  geom_errorbarh(aes(y="X", 
                 xmin = coefs_mod$coef - coefs_mod$sd*1.96, ## 95%
                 xmax = coefs_mod$coef +coefs_mod$sd*1.96), height=.02)+
  geom_vline(xintercept = 0)

Plot predictions with ggpredict

library(ggeffects)

plot(ggpredict(mod, "x"))

Replication exercice

In the folder replication, open the file Replication Kübler et al. (2025) - Exercice

Your task is to reproduce some of the analyses presented in the paper.

Group work

• Join in groups

• Create a Quarto document with your names, and title of your project

• Create the general structure of your report

• Try to import your data and do some first analyses

• Try to conduct a meaningful regression with your data for your project.